% Compute Results for Robustness Tables
%
clear all;
small = 1.0e-10;

% -- File Directories   
figdir = 'fig/';
outdir = 'out/';
matdir = 'mat/';
mcmcdir = 'mcmc/';

% Model Name -- Run once for each model
mod_name = '_baseline';
% mod_name = '_fgamma05';
% mod_name = '_fgamma15';
% mod_name = '_q_9';
% mod_name = '_q_23';
% mod_name = '_T68';
% mod_name = '_empdata';    % Y/L 

% % Compute Growth Rate Paths
if strcmp(mod_name,'_T68');
  n_t = 68;
elseif strcmp(mod_name,'_empdata');
  n_t = 68;
else;
  n_t = 118;
end;

% Load Posterior Elements
fstr = [mcmcdir 'sigma2_F_draws' mod_name];load(fstr);
fstr = [mcmcdir 'F_rho_draws' mod_name];load(fstr);
fstr = [mcmcdir 'stddev_del50U_draws' mod_name];load(fstr);
fstr = [mcmcdir 'halflife_draws' mod_name];load(fstr);
fstr = [matdir 'country_code']; load(fstr);
fstr = [mcmcdir 'g_matrix' mod_name];load(fstr);

sig_delta_a_draws = sqrt(sigma2_F_draws(2,:));
sig_m_draws = sqrt(sigma2_F_draws(1,:));
hl_f_draws = log(0.5)./log(F_rho_draws);
hl_pooled_draws = reshape(halflife_draws',[],1);
stddev_del50U_pooled_draws = reshape(stddev_del50U_draws',[],1);

pct = [1/6 0.5 5/6];
outfile_name = [outdir 'Table_4_Sensitivity' mod_name '.csv'];
fileID = fopen(outfile_name,'w');
fprintf(fileID,['Model: ' mod_name '\n\n']);
% Results for panel A
fprintf(fileID,'sigma_delta_a,');
tmp = 100*pctile(sig_delta_a_draws',pct)';
prtmat_comma(tmp,fileID,'%5.2f','\n');
fprintf(fileID,'sigma_m,');
tmp = 100*pctile(sig_m_draws',pct)';
prtmat_comma(tmp,fileID,'%5.2f','\n');
fprintf(fileID,'f-halflife,');
tmp = pctile(hl_f_draws',pct)';
prtmat_comma(tmp,fileID,'%5.2f','\n');

% Results for panel B
fprintf(fileID,'half-life,');
tmp = pctile(hl_pooled_draws,pct)';
prtmat_comma(tmp,fileID,'%5.0f','\n');
fprintf(fileID,'sigma_dc_50,');
tmp = pctile(stddev_del50U_pooled_draws,pct)';
prtmat_comma(tmp,fileID,'%5.2f','\n');

% % Panel C
fstr = [mcmcdir 'path_F_draws' mod_name];load(fstr);
fstr = [mcmcdir 'paths_U_draws' mod_name];load(fstr);
fstr = [matdir 'yp']; load(fstr);
fstr = [matdir 'pop']; load(fstr);
n_c = size(pop,2);
% % Rename;
u_draws = paths_U_draws;
f_draws = path_F_draws;
clear paths_U_draws;
clear path_F_draws;
% % 
% % Compute Growth Rate Paths
U_gr_paths_50 = squeeze(u_draws(n_t+50,:,:)-u_draws(n_t,:,:))/50;
U_gr_paths_100 = squeeze(u_draws(n_t+100,:,:)-u_draws(n_t,:,:))/100;
f_gr_paths_50 = (f_draws(n_t+50,:)-f_draws(n_t,:))/50;
f_gr_paths_100 = (f_draws(n_t+100,:)-f_draws(n_t,:))/100;
Y_gr_paths_50 = U_gr_paths_50 + repmat(f_gr_paths_50,n_c,1);
Y_gr_paths_100 = U_gr_paths_100 + repmat(f_gr_paths_100,n_c,1);
% 
fprintf(fileID,'100-yr F Growth,');
tmp = pctile(f_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f','\n');
% % 
% % Use population weights for growth rates
fprintf(fileID,'100-yr All growth,');
pop_g = pop(end,:)';     % All countries
wght = pop_g/sum(pop_g);
g_gr_paths_50 = wght'*Y_gr_paths_50;
g_gr_paths_100 = wght'*Y_gr_paths_100;
tmp = pctile(g_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
